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ABSTRACT 


The philosophy, history, operation, calibration of and some analyses with 
the ROAD (Rapid Orbit Analysis and Determination) program are described. 
This semi-numeric trajectory program integrates and analyses mean element 
variations for earth orbits with great efficiency. Through it's use, extensive 
zonal, resonant harmonic and earth tidal determinations have been made at 
Goddard Space Flight Center since 1969. 



pA 






iii 



CONTENTS 


Page 

INTRODUCTION 1 

HISTORY 3 

THE ROAD INTEGRATOR 4 

Geopotential Effects 4 

Third Body Gravity Effects (Direct) 9 

Drag Perturbations . 9 

Radiation Pressure Effects , 15 

Precession and Nutation . 17 

Solar and Lunar Tides . . 18 

Element Secular Rates and Accelerations 20 

ORBIT AND GEODETIC PARAMETER ESTIMATION 21 

RESULTS AND CALIBRATIONS 22 

SUMMARY 26 

REFERENCES 9.8 


pbj«edwg 








V 



THE ROAD PROGRAM 


INTRODUCTION 

At Goddard Space Flight Center, the most general computational system for 
analyzing long term gravitational and non-gravitational effects on satellite orbits 
has been the ROAD (Rapid Orbit Analysis and Determination) program. This 
program is a multiple arc/satellite orbit generator which can also estimate 
(from real data) a wide variety of the influential geodetic parameters. It uses 
mean Kepler elements as the usual data type. 

The ROAD orbit generator generally integrates numerically, orbit-averaged 
Kepler element variations determined from a number of sources. The variations 
due to the geopotential (to 40, 40), are determined from the right hand side of the 
Lagrange planetary equations in its standard form. Only those geopotential dis- 
turbances are usually integrated which have long term effects on the orbits, 
permitting time steps of the order of a day or more and extremely fast epheme- 
rides generation. Direct luni- solar gravity and the indirect luni-solar tidal ef- 
fects on the orbits are evaluated from a similar disturbing potential which 
selects only long period or orbit averaged terms for integration. Radiation 
pressure variations are evaluated from a quasi-potential for full sunlight con- 
ditions. When the orbit is partially in shadow, an or bit- aver aged force is evalu- 
ated by averaging the disturbance from shadow exit to shadow entrance. 



Precession and nutation effects of the movement of the earth’s polar axis are 
accounted for by integrating in a true of epoch inertial system. All forces are 

initially calculated in true of date coordinates and then rotated to true of epoch 

(using a precession-nutation rotation matrix). After integration, the true of 

epoch elements are then rotated back to true of date for comparison with real 

data which are usually given in these (latter) coordinates. 

Drag is the only perturbation not evaluated by a potential or quasi-potential. 
Therefore its element variations are not evaluated from the right hand side of 
the s tandar d planetary equations. Instead the gaussian form of the variation 
equations are employed with the normal, circumferential and radial components 
of the drag force evaluated from a one revolution J 2 perturbed trajectory through 
a model atmosphere. 

In its estimation mode, ROAD has the capability of solving for the following 
arc dependent parameters; The six initial Kepler elements, an offset or bias in 
the semi-major axis observations, up to 5 element rates and accelerations for 
each of the initial Kepler elements, and a drag and radiation pressure coefficient 
for each arc. In addition, it can solve for the following geodetic parameters, 
common to all the satellite arcs observed: The earth’s gaussian gravity constant 
and earth radius, geopotential harmonics to (40, 40), and earth tidal harmonics 
(love numbers) and their associated lag angles to 4th degree. ROAD determines 
these parameters by a Baysian least squares process using the Kepler element 
data as "observables." The partials of these observables with respect to the 
parameters are generally found by numerical integration of variation equations. 
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HISTORY 


The orbit generator of ROAD, and the basic idea of integrating mean or 
slowly varying Kepler elements, originated with the (resonant orbit) RESORB 
program developed by B. C. Douglas and G. S. Gedeon at TRW, Inc. in 1966 
[Gedeon, Douglas and Palmiter, 1967] . The original purpose of this program 
was to study the complex evolution of deeply resonant orbits for which no 
analytic theory exists. Later RE SORB, with only fixed geopotential and third 
body forces, was extended to include drag and became (in 1967) the satellite 
lifetime (Rapid Orbit Prediction) Program, ROPP [Wexler and Gedeon, 1967] , 
used extensively at Goddard Space Flight Center. Still later in 1968, ROPP 
was modified at Goddard Space Flight Center to give it a limited estimation 
capability (using secant partials) for state and geopotential recovery from long 
arcs of mean Kepler elements. In addition, the effects of radiation pressure 
were added, using Kaula's quasi-potential [Kaula, 1962], to account for signifi- 
cant changes in highly eccentric resonant orbits. 

The results of resonant geopotential determinations with this initial version 
of the ROAD (Rapid Orbit Analysis and Determination) program were published in 
Wagner, 1969, Wagner and Douglas, 1969 and Wagner, 1970b. In 1970, the long 
term rates due to the interaction of short period terms in the earth’s oblateness 
were added to the ROAD integrator so that a proper analysis of the evolution of 
close earth satellites could be made. Using this version of ROAD, preliminary 
results for zonal recovery from 2 satellites were published in Wagner, Putney and 
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Fisher , 1970. At the same time a new program was written using rigorous 
numerical integration of variation equations for partials and a complete Baysian 
(a priori) least squares scheme for parameter estimation [ Williamson , 1970, 
Guion , Lynn and Lynn, 1970 and Douglas, Dunn and Williamson, 1972] . Results 
from this new ROAD program were published in 1970 and 1972 on the determina- 
tion of both resonant and zonal geopotential terms from a large number of 
satellite orbits [ Wagner , 1970a, Wagner , 1972]. 

THE ROAD INTEGRATOR 

An Adams multistep integration scheme is used to solve the first order 
system of Lagrange's equations describing the satellite’s motion: 


d s. 

— = S. (s, p, t), s. (t = 0) = s. (0) Given ; i = 1, 2, . 


6 , 


(1) 


where the s are the satellite orbit's 6 state variables ; a, semi- major axis; e, 
eccentricity; I, inclination; w, argument of perigee; N, right ascension of the 
ascending node and M, the mean anomaly, and the p are geodetic parameters 
describing the force model. 

Geopotential Effects 

The rates S. for geopotential effects are derived from the standard form 
[Kaula, 1966] of (1): 

d a _ 2 3 R 

d t n a 3 M 
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where n = (m /a 3 ) 1/2 ; n being the satellite’s mean motion and ^ the earth’s 
gaussian gravity constant. R is the disturbing (non central) potential which, for 
the earth's gravitational field (fixed part as distinguished from time varying), 
has been expressed by Kaula [Kaula, 1966] in terms of the above Kepler ele- 
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is the earth's equatorial radius, 8 e is the Greenwich sidereal time, and 
(I) and G{, pq (e) are polynomials depending on inclination and eccentricity 
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alone, respectively. The and are the amplitudes and phase angles of the 
spherical harmonics of the potential, related to the ordinary C&S coefficients 
actually estimated in ROAD by 

Jt„ = * < c t * s L> ,/a 

= Tan -1 (Sp /C p )• 

The standard form of R e in terms of spherical harmonics is: 

03 l 't 

Re=j rZ] 2Z (t) p ^ sin ^ {c ^ m cosmX + s ^ sinmX> ’ <4) 

'{/ = 2 m= 0 

where r, 4> and X. are the satellite's distance from the center of mass of 
the earth, geocentric latitude and longitude; and the P| are the associated 
Legendre polynomials. 

The efficient functioning of the ROAD program depends upon the removal of 
all short periodic effects from the equations of motion. For the geopotential, 
this is accomplished by integrating only those terms (-C, m , p, q ) for which the 
argument of Sp is slowly varying. These long period terms fall into two 

'b mp q 

classes: 

For zonal harmonics (m - 0); % - 2p + q = 0, and for longitude dependent 
(resonant) harmonics (m / 0), £ - 2p + q = m/s are the indicial equations which 
must be satisfied. For the resonant case, s here is a rational number close to the 
actual mean motion (n) of the satellite in revolutions per day. The so called beat 
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period of the resonant orbit (with respect to the s commensurability) is the 
period it takes for the dominant argument of S to go through 360° when i - 2p + 
q = m/s. For near circular orbits this is generally the term for which q = 0. 

This scheme is equivalent to integrating mean, or orbit averaged Kepler 
elements and results in the proper interaction of all long-periodic effects with 
each other. However, interactions of short-periodic terms (-1 - 2p + q / m/s, for 
m / 0 and •£ - 2p + q / 0, for m = 0) with each other and with long-periodic 
terms are lost. But the only case where this neglect is significant are the long 
periodic effects caused by the mutual interaction of the short periodic effects 
due to the earth's oblateness (J 2 0 ). When only long period terms are included, 
the ROAD integrator adds the following element rates derived from the second 
order oblateness potential due to Brouwer [Brouwer, 1959]; 
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(5) 


where the Delaunay variables are: L = (jj. a) 1/2 , G = L (1 - e 2 ) 1/7 and 
H = G cos I; and n Q = m 2 /L 3 , y 2 - k 2 /L 4 , k 2 = J 2 (a e ) 2 /2 t and J 2 = - J 2 ,o* 
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Third Body,, Gravity Effects (Direct) ;./ v. 

For the 3rd body perturbations, an entirely analogous development has been 
made in Kaula 1962, and is used in ROAD. In tUs case, the perturbations, can be 
, sorted out according to frequency in much the same manner as the geopotential 
perturbations. The long periodic terms of the 3rd body disturbing function are 
the only ones presently coded in ROAD and have the form - .. 

= f (a, e, I, a*, e*, I*) cos [(n - 2 h) w* - (n - 2 p) w + (n - 2 h + j) M* 

+ m (N - N*)], • 

where n here is the degree of the third body potential (2 £ n £ 4, as coded, and 
where starred quantities refer to orbit elements of the disturbing body. Analogous 
to the geopotential case; 0£h£n, 0£p£n and 0 <; min. Since f is proportional 
to (e*) I j I, only - 4 i j £ 4 is coded and can be limited further on option. A 
further option is to ignore the less than monthly lunar terms where n - 2h + j ^ 0. 

Drag Perturbations 

The method of computation of the long periodic and secular variations of 
the elements due to atmospheric drag has undergone extensive development in 
ROAD. Originally, analytic expressions similar to those given in King-Hele, 

1964 were used, but when these expressions were modified to include non- 
linearity of scale height and time dependence, the complexity became very great, 
and separate developments for low and high eccentricities were required 
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[Wexler and Gedeon, 1967] . Starting in 1971 this method was abandoned in favor 
of a simple, purely numerical scheme (B. Chovitz, personal communication) 
originally including only the effects on a and e. 

The integrals that give the average time rates -of- change of the elements due 
to drag are evaluated numerically by a 9-point Gauss quadrature integration. At 
each integration point it is necessary only to compute the altitude of the satellite, 
and evaluate the density of the atmosphere from any atmosphere model. The 
ROAD program currently uses either the Jacchia 1968 or 1971 atmosphere 
models [Diamante and Der, 1972]. 

The element rates are computed using the Gaussian form of the equations 
of motion, which is expressed in terms of radial <F ), transverse (F t ), and 
normal forces (F ). The drag force is the usual simple model: 


where 



v 

r 


C Q is the satellite drag coefficient 

A is the projected cross-sectional area of the satellite normal to v 

r 

m is the mass of the satellite 

P is the density of the atmosphere at the satellite position, and 
v r is the velocity vector of the satellite relative to the atmosphere. 

The most significant problem in computing the correct drag is the 
calculation of an accurate spheroid height of the satellite. This is done on the 
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osculating J 2 perturbed orbit over which the drag is averaged. The formula used 
was developed in Ko 2 ai, 1959 and modified to use the mean semi- major axis 
instead of his "geometric" semi-major axis. That is: 


A r = 




1 + — (l - / I - e- ) cos f 
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— - sin f sin 2 w 

1 - e 2 J 


where p « a(l - e 2 ) and f is the true anomaly. 


Thus the satellite's radius vector has magnitude r + Ar where r is computed from 
the two- body formula. 

The spheroid height is then obtained by subtracting the spheroid radius of 
the Earth at the latitude under consideration: 


S E .,«h = a . [l- ( ft f f2 ) sin2 *' +yf 2 >in‘*'l : 

where in this particular formula, 

f is the flattening of the Earth, and 
4 >‘ is the satellite's geocentric latitude. 
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The atmosphere is assumed to rotate with the Earth. The rotation rate of 
the Earth resolves into a component ( 9 e cos I) normal to the orbit and a component 
(6 e sin I) perpendicular to the angular momentum vector, h, of the satellite, which 
is perpendicular to the vector, N, pointing toward the ascending node. 

The satellite's velocity vector is given by 

T = (r, r f, 0) 

in radial, transverse (or circumferential), and normal components, where r is 
the satellite's position vector. 

The relative velocity vector is then 

v r - [ r , rf-(9 e r cos I, 6 r sin I cos (w + f)] . 

From the vis viva or energy integral, 

7 ■ 7 = r 2 + r 2 f 2 = ft - 

Also, the angular momentum h is given by 

h = r 2 f = /jU. a (1 — e 2 ) . 

Hence: 



v r = M 


(2 1 \ 2 h 4 _ 

1 Q cos I 

\ r a/ r e 


+ r2 [cos 2 I + sin 2 I cos 2 (w + f)] . 
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This last formula is used to compute the magnitude' of the relative -velocity vector. 

The Gaussian form of the variation equations for the perturbations of the 
Kepler elements are: 
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For the drag case- being considered: 


F r = - F d r 






F ,=-F d (rf - 0 e r co S I)" *5 


F n = - F d (8 e r sin I cos (w + f) 
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Using the usual two body formulae for relating Kepler variables, the element 
rates are simplified to: 


da ^ 2c , /l + e 2 + 2 e cos f p c os ^ 
dT = ~ 2a \ p h 


d e 
d t 
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These equations are now averaged over one revolution in True Anomaly to ob- 
tain the mean element rates which are added (on option) to the other rates used 
in the ROAD integrator. This averaging is done numerically by a 9 point Gauss 
quadrature method with the middle point placed at perigee to insure an accurate 
computation of the maximum drag. 
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Radiation Pressure Effects 


The scheme used to calculate the long-periodic radiation pressure effects 
is that given in Kaula, 1962. He notes that for a satellite entirely in sunlight, a 
quasi-potential (for the orbit averaged disturbance) can be developed for radia- 
tion pressure. Kaula’s quasi-potential is used as is in the Lagrange Planetary 
equations for sunlight conditions. However, if a shadow is present, an integra- 
tion of the derivatives of the full potential from exit point to entry point is 
made to obtain the orbit-averaged disturbance. The values of the eccentric 
anomaly at these points are obtained by solution of the quartic equation given in 
Kaula’s paper. The satellite is assumed to be spherical, and the earth’s shadow 
cylindrical. The radiation force f r (actually, the acceleration) is given as (I/C) 
(A/M) C R where I is the Solar Flux in space (1.37 x 10 6 erg/cm 2 -sec) and C is 
the velocity of light. The original subroutine for these effects in ROAD was 
supplied by Bernard Chovitz [B. Chovitz, Personnal Communication, 1970] i 

For a shadowless orbit, the long-term (orbit averaged) potential is as given 
in Kaula, 1962: 

R = - f r (3 a e/2) [cos 2 (1/2) sin 2 (e/2) cos (w + N + X Q ) - 
+ cos 2 (1/2) cos 2 (e/2) cos (w + N - \ @ ) 

+ sin 2 (1/2) cos 2 (e/2) cos (w - N + \ ) 

- 1/2 sin I sine cos (w + \ 0 ) + 1/2 sin I sin e cos (w - \ 0 )] •, 

where e and \ are the inclination of the ecliptic ( 23.45°) and the sun's true 
longitude. 
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The derivatives of the quasi-potential with respect to the orbit elements are 
then used in the standard form of the Lagrange equations, (2), to form the 
element rates from the radiation pressure for the shadowless orbit. 

For the shadowed orbit, the long term, orbit averaged, derivatives are given 
in terms of the exit and entry eccentric anomalies, E 0 and E 1 , as: 
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where R gq is a rotation matrix relating the position of the satellite in an orbital 
plane coordinate system to its position with respect to a coordinate system point- 
ing at the sun. \R is defined as the successive rotations: 

0 — s q 


where, 


R 3 (0*B 3 C-N)-r 1 (- I)-R 3 (- w), 



Z cos x s in x 0 
— j sin x cos x 0 

^ 0 0 1 

R n and R 12 are elements of R sq . 

Precession and Nutation 
The effects of precession and nutation of the earth's polar axis are ac- 
counted for in ROAD by maintaining the integration in an inertial system with 
respect to the equinox and true equator at epoch (TE). All the forces are initially 
calculated in a system referred to the equinox and true equator of date (TD), 
corresponding to the observed data. The forces (actually the Kepler rates) are 
then rotated (by a precession-nutation matrix) to refer to the TE system where 
the integration proceeds. The integrated state is rotated back to the TD system 
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at the end of each Integration step for calculation of new rates. The chain of 
rotations actually used to obtain the rates in the TE system is: 

^ S TE / 'd S TE ^ X TE ^ ’S’D \ ^ S TD 

dt \**t* B *rn dt 

where x is the cartesian position-velocity vector corresponding to the Kepler 
elements s in the indicated coordinate system. The matrix dx/d s is documented 
on page 68 of Kaula's book [Kaula, 1966]. The precession-nutation matrix 
cix /Bx is described in volume I of the GEODYN Documentation, (Sept. 30, 
1972, Goddard Space Flight Center NAS-5-11735 Mod. 65, PCN 550W-72416) and 
provides for Simon Newcomb's description of the precession and Edgar Woolard's 
nutation. 

Solar and Lunar Tides 

The ROAD program treats the tidal perturbations in a manner analogous to 
all other gravitational perturbations. The potential is expanded in terms of 
Kepler elements, and the long periodic terms only are selected for differentiation 
in the computation of the Kepler element rates from the standard form of the 
Lagrange equations (2). The form of the potential is due to Kaula t Kaula, 1969], 
and includes provision for latitude dependent tidal amplitudes (Love numbers) 
and phase lags (of the potential bulge from the 3rd body- earth line). 

The potential that is used in ROAD is equation (22) in Kaula, 1969: 
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The subscripted k's are the love number's and the <='s are the phase lags of the 
tidal potential. In the above potential, no short period effects (in the satellite's 
mean anomaly) are considered. Both K and Q „ terms (the dominant ones 
for a distant 3rd body; l = 2) are defined in tables II and I of Kaula's 1969 paper. 
The F and G functions are the usual ones for satellite orbits defined in Kaula's 
book [Kaula, 1966 ). The latitude dependence of the love number k ^ is defined by 

k e = L k ^H P H0 C sin ^)> 

h - 

and the phase lag by 


m P <0 P nn ( sin <£) 


where thep's are the ordinary Legendre polynomials. The quantity (fye^is 
the sum 
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The longitude arguments are: 



e Q 

n ^ 


knsm' 


V, - v? = [(k - 2 j) w + m N] - H - 2 p) w* + {i - 2 p + q)M* + m N*], 

kmj ( 2) - k) 'Vmpq 

where the starred elements belong to the 3rd body. 

Element Secular Rates and Accelerations 

The ROAD program also can include the effect of up to 5 arbitrary element 
rates and accelerations (derivatives) of the six initial Kepler elements El of 
each arc. Thus, the elements E at time t are expanded in a Taylor series in the 
arbitrary derivatives E^ at the epoch time t 0 : 


n=0 


(t _ t) n r l 

“ E U (t„) 


In addition, the program can propagate a constant bias in each arc of semi-major 
axis observations, to satisfy the equation: 


a (OBSERVED) = a (CALCULATED) + a (BIAS). 

This bias is often needed to account for the different definitions of the mean 
semi-major axis between ROAD and the observations used. 
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ORBIT AND GEODETIC PARAMETER ESTIMATION 

The logic of ROAD in its estimation (or differential correction) inode follows 
the (now universally used) scheme of partitioned normals [Kaula, 1966]* Param- 
eters are separated into "arc" parameters; that is, parameters unique to each 
orbital arc, aiid common parameters, or parameters common to all orbit arcs. 

In a multi-arc solution with common parameters, the common parameter matrix 
is first solved and then the arc matrices are successively solved by back- 
substitution. 

Of course parameter and orbit estimation requires partial derivatives. In 
ROAD these are obtained by numerically integrating variational equations. The 
equations of motion have the form 

4*r = € f (a, e, i, N, w, M, t) 
a t 

where s is an array Kepler elements, and e is a parameter to be estimated. Thus 
formation of the variational equations is elementary. The rigorous variational 
equation for a parameter € that is numerically integrated is: 



f (a,, e, i, N, w, M, t) 
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However, because of the small interaction of effects, the second term on the right 
hand side is ignored for radiation pressure, Earth tides, and even drag parameter 
estimation. 

In parameter estimation the ROAD program utilizes a Baysian least 
squares scheme; i.e., a-priori information is required on all adjusted param- 
eters. However the program accepts only a-priori standard deviations. No 
covariances are allowed (they are assumed to be zero). In ROAD, a-priori 
information is used to control the conditioning of a problem and also the amount 
by which parameters are allowed to adjust. But because the full covariances 
information is not given, this adjustment control is not always rigorous. 

RESULTS AND CALIBRATIONS 

The ROAD program has undergone extensive calibration on long arcs of 
simulated osculating element data in the course of its development [i.e., Wagner , 
Putney and Fisher, 1970 and Gulon, Lynn and Lynn, 1970 1 While these 
calibration results have generally been successful they point up the need 
for accurate osculating to mean element conversion when ROAD operates in its 
rapid mean element mode. Since the most accurate observed elements on real 
orbits are osculating (i.e.; derived by numerical integration from complete force 
models), the availability of a good converter is not just academic. The results 
presented here (on both real and simulated data) were obtained with the use of a 
combined analytic-numerical converter discussed in Douglas, Marsh and Mullins, 
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1972. Previous ROAD preprocessing of osculating elements [Wagner, 1972] 
removed short period effects using the analytic theory of Brouwer , 1959. The 
new method first removes analytically all short period effects due to the geo- 
potential through (4, 4) tKaula, 1966], The remaining effects are averaged 
numerically over an orbit revolution with respect to a secularly precessing 
ellipse. 

The satellites used for the real data calibration were GEOS-1, GEOS-2, 
and PEOLE. GEOS 1 has a near critical inclination, and so experiences rather 
large perturbations from odd zonal harmonics. The orbit is also slightly eccen- 
tric (e = 0.07) so that radiation pressure perturbations are important. Drag, 
although very small, is detectible. 

GEOS-2 has a more nearly circular orbit than GEOS-1, and is thus rela- 
tively much less affected by radiation pressure. However, the high (106°) in- 
clination of the orbit causes a slow node rate (0^6 /day) and resulting large tidal 
perturbations exceeding 10 arc seconds in the inclination and 30 arc seconds in 
the node. 

PEOLE is a low inclination (15°), low perigee satellite with rapid node and 
perigee rates. Drag is very high, the semi- major axis decaying about 30 m per 
day. All of these satellites are significantly perturbed by the sun, moon, and 
the effects of precession and nutation. 

The orbit of PEOPLE is a good test of the ability of ROAD to represent 
properly the effects of atmospheric drag. Figure 1 shows the observed and 
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computed semi-major axes of PEOLE obtained from a ROAD orbit determination. 
The time covered is the year 1971. The mean elements were determined from 4 
day orbital arcs of minitrack and laser data as part of the International Satellite 
Geodesy Experiment (ISAGEX) [Marsh, Douglas and Klosko, 1971]. Note the 
close agreement of observed and computed elements using the Jacehia 1968 
model atmosphere. Subtle variations are visible and are properly modeled.* 

The overall error for the year is only about 5%, and appears to be long-periodic. 
It is suspected that the Jacchia atmosphere of 1968 may have an error in the 
semi-annual variation of density, so a small unmodeled variation is not 
surprising (see also Wagner, 1972). 

Figure 2 shows the mean eccentricity of PEOLE, also obtained from the 
ISAGEX data. The periodic variation due to the odd zonal harmonics is clearly 
visible, as is the secular decrease due to drag. Figure 3 shows the residual 
eccentricity with drag modeled as above and zonal harmonics from the SAO 
1969 Standard Earth [Gaposchkin and Lambeck, 1970], Note that drag effects 
are removed essentially perfectly (the residuals have no clear secular trends) 
but a residual odd zonal harmonic effect remains. Since the SAO 1969 zonal 
harmonics included no low inclination satellites, such a residual effect is to be 
anticipated. To verify that this is not due to error in ROAD, the eccentricity of 
PEOLE was obtained from a complete "Cowell "-type, numerical integration 
of the geopotential to (4, 4). Mean elements were then prepared from 

* Where observed and computed points are coincident, only the observed is printed. 
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this trajectory by the method in Douglas , Marsh, and Mullins (1972), Figure 4 
shows the variation in eccentricity produced by C, n (predominantly) and Figure 5 
shows the ROAD fit to this data using the same value of C 3 The residuals ap- 
pear random and have an rms of only 0.1% of the amplitude of the eccentricity 
variation shown in Figure 4. 

GEOS-1 offers a good opportunity to study the ROAD radiation pressure 
model; Drag is small for GEOS-1 (but not negligible). Resonance produces a 
small, significant effect, but with a period of less than a week. Because of the 
substantial eccentricity of the orbit, radiation pressure causes the dominant 
perturbation of the semi- major axis during orbit shadow periods. Fluctuations 
of 5-15 meters over several months are typical. These can be studied in detail 
because very accurate mean elements for GEOS-1 are available. The mean 
elements obtained in Douglas, Marsh, and Mullins (1972) from 2 day orbital 
arcs of optical data have a precision of 25 cm in the semi-major axis (see 
figure 6). 

The decay of semi-major axis of GEOS 1, due to drag, amounts to 5-10 m per 
year and is most readily seen during a no-shadow period, such as July-August 
1966. Because of these significant drag effects, it was found necessary to ac- 
count for both drag and radiation pressure perturbations on GEOS-1 . 

The radiation pressure model in ROAD does not include the Earth Albedo 
effect. Thus a small error in the modeling of radiation pressure is possible. 

The period March 11 - May 15, 1966 on GEOS-1 has been studied intensively for 
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determination of tidal parameters where, of course, accurate modeling of radia- 
tion pressure is essential. During this period the semi-major axis is perturbed 
in a non-linear fashion by radiation pressure by about 5-6 meters, and in a linear 
way from drag by about 60 cm. Figure 7 shows the residuals obtained in the mean 
semi-major axis for this period by ROAD when C D and C R are adjusted to the 
values C = 3,1, and C„ = 1.52. The rms of the fit is less than 25 cm. 

Another important test of ROAD concerns Earth tides. Figure 8 shows the 
residuals in the inclination of GEOS 2 over a 600 day period with respect to mean 
elements obtained from 2 day arcs of optical tracking data. Earth tides were not 
modeled in obtaining these residuals. Note a long period unmodeled perturbation 
of amplitude 5" in the inclination. Figure 9 shows the results with the Love 
number k 2 =0.30 modeled in ROAD, The residuals are almost random. This 
value of k 2 also substantially explains a similar, substantial, solar tidal effect 
(of the same long period) in the node of GEOS 2. This value of the Love number 
is in reasonable agreement with values derived from other orbits which range 
from 0.25 to 0.35. [ Kaula, 1969] . 

SUMMARY 

The ROAD program can efficiently analyze long arcs of mean elements for a 
wide variety of geodetic effects. Since 1969, the program has been used extensively 
to determine and evaluate both resonant and zonal geopotential harmonics from 
a large number of orbits. It has also been used to investigate possible variations 
in low order earth-gravity ''constants” [Wagner, 1973], 
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Currently, the program is also being used at Goddard Space Flight Center to 
study the time varying geopotential due to sun and moon induced earth tides, 
through their effects on a large number of well tracked satellite orbits 
[Douglas and Marsh, 1973] . 
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Figure 1. The Semi-major Axes of PEOLE Obtained From a Road Orbit Determination. 
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Figure 2. The Eccentricity of PEOLE. 
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Figure 3. Eccentricity Residuals for PEOLE. 
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Figure 4. The Mean Eccentricity From a Simulated PEOLE Trajectory. 
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Figure 5. Residuals in Eccentricity Computed by ROAD From a Simulated PEOLE Trajectory. 
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Figure 7. Residuals in the Semi-Major Axis from a ROAD Orbit 
Determination of GEOS-1 Adjusting C and CL . 
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Figure 8. Residuals in the Inclination of GEOS-2 from a ROAD Orbit Determination 
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Figure 9. Residuals in the Inclination of GEOS-2 Modeling Tides at K 2 = .30. 




